##################################################
## Replication materials for 
## "Competent legislators or mere pawns? 
## Experimental evidence of attitudes toward gender quota politicians"
## by Carolyn Barnett, Alexandra Blackman, and Marwa Shalaby

## 6. Robustness Checks 
## This file creates Appendix Figures A.7 to A.11 and Appendix Tables A.5 to A.12

## **NOTE: Run file "01_setup.R" first to load necessary packages and data**

# MISSINGNESS ---------------

## Data quoted in text -------------

# absolute # of NA responses
table(is.na(morvign$vignette_competence)) # 21 NA
21/927 # 2.3%

table(is.na(morvign$vignette_pawnlike)) # 34 NA
34/927 # 3.7%

table(is.na(morvign$vignette_cooperation)) # 25 NA 
25/927 # 2.7%

# identify people who ever non-respond
morvign <- morvign %>%
  mutate(
    comp_na = ifelse(is.na(vignette_competence), 1, 0),
    coop_na = ifelse(is.na(vignette_cooperation), 1, 0),
    pawn_na = ifelse(is.na(vignette_pawnlike), 1, 0)
  )

table(morvign$comp_na)

na_data <- morvign %>%
  dplyr::filter(comp_na == 1 | 
                  coop_na == 1 |
                  pawn_na == 1)
length(unique(na_data$QID)) # 47 unique respondents
47/927 # 5.1 %

## FIGURE A.7: Correlates of Non-Response in Vignette Outcomes ----------

# create indicator in dataset for ever non-responding
na_resps <- unique(na_data$QID)
morvign <- morvign %>%
  dplyr::mutate(
    vig_nas = ifelse(QID %in% na_resps, 1, 0)
  )
table(morvign$vig_nas)

# create education variable for regression
morvign <- morvign %>%
  dplyr::mutate(education = case_when(
    edu_no_formal == 1 ~ 1,
    edu_elem_prep == 1 ~ 2,
    edu_second == 1 ~ 3,
    edu_above_second == 1 ~ 4
  ))

# Regression with everything and plot for appendix
missing <- lm(vig_nas ~ female + age + education +
                ses_subj + religiosity + govt_satis +
                voted_rni + voted_pjd + voted_usfp +
                voted_istiqlal,
              data = morvign
)

cm <- c('female' = 'Female',
        'age' = 'Age', 
        'education' = 'Education',
        'ses_subj' = 'Socioeconomic Status',
        'religiosity' = 'Religiosity',
        'govt_satis' = 'Govt Satisfaction',
        'voted_rni' = 'Voted RNI',
        'voted_pjd' = 'Voted PJD',
        'voted_usfp' = 'Voted USFP',
        'voted_istiqlal' = 'Voted Istiqlal')

modelsummary::modelplot(missing,
                        coef_omit = 'Interc',
                        coef_map = cm) +
  geom_vline(xintercept = 0, linetype = 3)

ggsave("missingness_corr_vignette.png", plot = last_plot(),
       path = "Figures/", limitsize = F,
       bg="white",
       width = 5.4, height = 5.7, units = "in", dpi = 300)


# MANIPULATION CHECK SAMPLE --------------------
## TABLE A.5: Manipulation Check Responses by Treatment Group -----------

# calculate frequencies 
manip_results <- data.frame(table(morvign$vignette_election, 
                                  morvign$manip_correct)) %>%
  mutate(Response = case_when(
    Var2 == 0 ~ "Incorrect",
    Var2 == 1 ~ "Correct",
    Var2 == 9 ~ "Don't Know or Refused"
  )) %>%
  dplyr::select(-Var2) %>%
  pivot_wider(names_from = "Response",
              values_from = "Freq") %>%
  relocate(Correct, .before = "Incorrect")

# calculate percentages
manip_results_prop <- data.frame(prop.table(table(morvign$vignette_election, 
                                                  morvign$manip_correct),
                                            margin = 1)) %>%
  mutate(Response = case_when(
    Var2 == 0 ~ "Incorrect",
    Var2 == 1 ~ "Correct",
    Var2 == 9 ~ "Don't Know or Refused"
  )) %>%
  dplyr::select(-Var2) %>%
  pivot_wider(names_from = "Response",
              values_from = "Freq") %>%
  relocate(Correct, .before = "Incorrect")

# combine results for table
manip_print <- rbind(manip_results[1,],
                     manip_results_prop[1,],
                     manip_results[2,],
                     manip_results_prop[2,],
                     manip_results[3,],
                     manip_results_prop[3,])

manip_print <- round(manip_print[,2:4], digits = 2)

print(xtable(manip_print))


## FIGURE A.8: Effect of politican gender and mode of election: manipulation check sample ---------

# Restrict sample to those who answered correctly
mv_manip <- morvign %>%
  filter(manip_correct == 1) # 662 obs

# Transform to long - NOTE: replace "morvign" in original code with "mv_manip" then run as written
morvign_long <- gather(mv_manip, measure, value, 
                       c(vignette_competence:vignette_pawnlike), factor_key=TRUE)

# Recode variables
morvign_long$vignette_election <- factor(morvign_long$vignette_election)
morvign_long$vignette_mode <- factor(morvign_long$vignette_mode)
morvign_long$vignette_work <- factor(morvign_long$vignette_work)
morvign_long$measure <- factor(morvign_long$measure, levels= c("vignette_competence", 
                                                               "vignette_cooperation", 
                                                               "vignette_pawnlike"))

vign_election_means <- summarySE(morvign_long, measurevar="value", 
                                 groupvars=c("measure","vignette_election"), na.rm= T)
vign_election_means

vign_election_main <- vign_election_means %>% 
  filter(!is.na(vignette_election)) %>%
  filter(measure != "vignette_cooperation")

##### CREATE MAIN OUTCOMES GRAPH (FIGURE 1)

ggplot(vign_election_main, aes(x=measure, y=value, fill=vignette_election)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=value-ci, ymax=value+ci),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9),  col = "black")+
  ylab("Outcome mean") + 
  scale_fill_manual(values = c("#666666", "#999999", "#CCCCCC")) + 
  coord_cartesian(ylim = c(3.5, 6.5)) +
  scale_x_discrete(labels=c("Competence","Pawnlike")) +
  labs(x = "Outcome measures", 
       fill = "Politician Gender &\nMode of Election", 
       caption = "")

ggsave("vign_election_main_manip.png", plot = last_plot(),
       path = "Figures/", limitsize = F,
       width = 7.5, height = 4.5, units = "in", dpi =300)


## TABLE A.6: Differences in Means for Type of Collaboration: Manipulation Check Sample --------

###### Create just quota women dataset
morvign_wq <- mv_manip %>% 
  filter(vignette_election=="Woman/quota")

# create tidy table with t-tests for variables of interest
morvign_women_work_t_tests <- bind_rows(tidy(t.test(vignette_competence ~ vignette_work, morvign_wq)),
                                        tidy(t.test(vignette_cooperation ~ vignette_work, morvign_wq)),
                                        tidy(t.test(vignette_pawnlike ~ vignette_work, morvign_wq))
)

# prepare table for latex estimate1 = general women & estimate2 = quota women
morvign_women_work_t_tests <- morvign_women_work_t_tests %>%
  add_column(variable = c("Competence",
                          "Cooperation",
                          "Pawnlike")) %>%
  dplyr::select(variable,estimate1,estimate2,estimate,p.value) %>%
  mutate_at(vars(estimate1, estimate2, estimate), funs(round(., 3))) %>%
  mutate_at(vars(p.value), funs(round(., 3))) %>%
  mutate(estimate = -estimate) %>%
  dplyr::rename("Outcome" = variable,
                "Mean (With his/her party)" = estimate1,
                "Mean (With another party)" = estimate2,
                "Difference" = estimate,
                "p-value" = p.value)

stargazer(morvign_women_work_t_tests,
          type="latex",
          summary = F,
          rownames = F,
          single.row = T)

## FIGURE A.9: Mean Response by Collaboration Treatment Among Quota Women: Manipulation check sample -----------

# *Note: need to create "morvign_long" from "mv_manip" in the section above for Figure A.8 first

# dataset of only quota women
vign_quota_collab <- morvign_long %>% 
  filter(!is.na(vignette_election)) %>%
  filter(vignette_election == "Woman/quota")

vign_quota_collab_means <- summarySE(vign_quota_collab, 
                                     measurevar="value", 
                                     groupvars=c("measure","vignette_work"), na.rm= T)
vign_quota_collab_means

vign_quota_collab_means <- vign_quota_collab_means %>% 
  filter(!is.na(vignette_work)) %>%
  filter(measure != "vignette_cooperation")


##### CREATE GRAPH OF COLLABORATION TYPE (QUOTA WOMEN ONLY)

ggplot(vign_quota_collab_means, aes(x=measure, y=value, fill=vignette_work)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=value-ci, ymax=value+ci),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9),  col = "grey70")+
  ylab("Outcome mean") + 
  scale_fill_manual(values = c("#999999", "#333333")) + 
  coord_cartesian(ylim = c(4.1, 6.5)) +
  scale_x_discrete(labels=c("Competence","Pawnlike")) +
  labs(x = "Outcome measures", fill = "Type of Collaboration", 
       caption = "")

ggsave("vign_quota_collab_means_manip.png", plot = last_plot(),
       path = "Figures/", limitsize = F,
       width = 7.5, height = 4.5, units = "in", dpi =300)

## TABLE A.7: Effect of Working with Another Party on Competence Rating: Manipulation Check Sample, All Treatment Arms --------

# Effect of type of collaboration on competence
vign_h2_comp1 <- lm(vignette_competence ~ vignette_work, 
                    data = mv_manip[mv_manip$vignette_election=="Man/general",])

vign_h2_comp2 <- lm(vignette_competence ~ vignette_work, 
                    data = mv_manip[mv_manip$vignette_election=="Woman/general",])

vign_h2_comp3 <- lm(vignette_competence ~ vignette_work, 
                    data = mv_manip[mv_manip$vignette_election=="Woman/quota",])

# Print table
stargazer(vign_h2_comp1,vign_h2_comp2,vign_h2_comp3,
          title = "Effect of Working With Another Party on Competence Rating: Manipulation Check Sample, All Treatment Arms",
          font.size = "scriptsize",
          keep.stat = c("n","adj.rsq"),
          model.names = F,
          multicolumn = TRUE)


## TABLE A.8: Effect of Working With Another Party on Pawn-Like Rating: Manipulation Check Sample, All Treatment Arms ---------

# Effect of type of collaboration on pawnlike
vign_h2_pawn1 <- lm(vignette_pawnlike ~ vignette_work, 
                    data = mv_manip[mv_manip$vignette_election=="Man/general",])
vign_h2_pawn2 <- lm(vignette_pawnlike ~ vignette_work, 
                    data = mv_manip[mv_manip$vignette_election=="Woman/general",])
vign_h2_pawn3 <- lm(vignette_pawnlike ~ vignette_work, 
                    data = mv_manip[mv_manip$vignette_election=="Woman/quota",])

# Output table format
stargazer(vign_h2_pawn1,vign_h2_pawn2,vign_h2_pawn3,
          title = "Effect of Working With Another Party on Pawn-Like Rating: Manipulation Check Sample, All Treatment Arms",
          font.size = "scriptsize",
          keep.stat = c("n","adj.rsq"),
          model.names = F,
          multicolumn = TRUE)



# NON-STRAIGHLINERS SAMPLE ------------------

## FIGURE A.10: Effect of politician gender and mode of election: Non-Straightliners Sample -------

# Restrict sample to those who did not straightline
mv_straight <- morvign %>%
  dplyr::filter(straightliner == 0) # 879 respondents

# Transform to long - NOTE: replace "morvign" in original code with "mv_straight" then run as written
morvign_long <- gather(mv_straight, measure, value, 
                       c(vignette_competence:vignette_pawnlike), factor_key=TRUE)

# Recode variables
morvign_long$vignette_election <- factor(morvign_long$vignette_election)
morvign_long$vignette_mode <- factor(morvign_long$vignette_mode)
morvign_long$vignette_work <- factor(morvign_long$vignette_work)
morvign_long$measure <- factor(morvign_long$measure, levels= c("vignette_competence", 
                                                               "vignette_cooperation", 
                                                               "vignette_pawnlike"))

vign_election_means <- summarySE(morvign_long, measurevar="value", 
                                 groupvars=c("measure","vignette_election"), na.rm= T)
vign_election_means

vign_election_main <- vign_election_means %>% 
  filter(!is.na(vignette_election)) %>%
  filter(measure != "vignette_cooperation")

##### CREATE MAIN OUTCOMES GRAPH (FIGURE 1)

ggplot(vign_election_main, aes(x=measure, y=value, fill=vignette_election)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=value-ci, ymax=value+ci),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9),  col = "black")+
  ylab("Outcome mean") + 
  scale_fill_manual(values = c("#666666", "#999999", "#CCCCCC")) + 
  coord_cartesian(ylim = c(3.5, 6.5)) +
  scale_x_discrete(labels=c("Competence","Pawnlike")) +
  labs(x = "Outcome measures", 
       fill = "Politician Gender &\nMode of Election", 
       caption = "")

ggsave("vign_election_main_straight.png", plot = last_plot(),
       path = "Figures/", limitsize = F,
       width = 7.5, height = 4.5, units = "in", dpi =300)


## TABLE A.9: Difference in Means for Type of Collaboration: Non-Straightliners Sample --------------

###### Create just quota women dataset
morvign_wq <- mv_straight %>% 
  filter(vignette_election=="Woman/quota")

# create tidy table with t-tests for variables of interest
morvign_women_work_t_tests <- bind_rows(tidy(t.test(vignette_competence ~ vignette_work, morvign_wq)),
                                        tidy(t.test(vignette_cooperation ~ vignette_work, morvign_wq)),
                                        tidy(t.test(vignette_pawnlike ~ vignette_work, morvign_wq))
)

# prepare table for latex estimate1 = general women & estimate2 = quota women
morvign_women_work_t_tests <- morvign_women_work_t_tests %>%
  add_column(variable = c("Competence",
                          "Cooperation",
                          "Pawnlike")) %>%
  dplyr::select(variable,estimate1,estimate2,estimate,p.value) %>%
  mutate_at(vars(estimate1, estimate2, estimate), funs(round(., 3))) %>%
  mutate_at(vars(p.value), funs(round(., 3))) %>%
  mutate(estimate = -estimate) %>%
  dplyr::rename("Outcome" = variable,
                "Mean (With his/her party)" = estimate1,
                "Mean (With another party)" = estimate2,
                "Difference" = estimate,
                "p-value" = p.value)

stargazer(morvign_women_work_t_tests,
          type="latex",
          summary = F,
          rownames = F,
          single.row = T)

## FIGURE A.11: Mean Responses by Collaboration Treatment Among Quota Women: Non-Straightliners Sample --------

# *Note: need to create "morvign_long" from "mv_straight" in the section above for Figure A.10 first

# dataset of only quota women
vign_quota_collab <- morvign_long %>% 
  filter(!is.na(vignette_election)) %>%
  filter(vignette_election == "Woman/quota")

vign_quota_collab_means <- summarySE(vign_quota_collab, 
                                     measurevar="value", 
                                     groupvars=c("measure","vignette_work"), na.rm= T)
vign_quota_collab_means

vign_quota_collab_means <- vign_quota_collab_means %>% 
  filter(!is.na(vignette_work)) %>%
  filter(measure != "vignette_cooperation")


##### CREATE GRAPH OF COLLABORATION TYPE (QUOTA WOMEN ONLY)

ggplot(vign_quota_collab_means, aes(x=measure, y=value, fill=vignette_work)) + 
  geom_bar(position=position_dodge(), stat="identity") +
  geom_errorbar(aes(ymin=value-ci, ymax=value+ci),
                width=.2,                    # Width of the error bars
                position=position_dodge(.9),  col = "grey70")+
  ylab("Outcome mean") + 
  scale_fill_manual(values = c("#999999", "#333333")) + 
  coord_cartesian(ylim = c(4.1, 6.5)) +
  scale_x_discrete(labels=c("Competence","Pawnlike")) +
  labs(x = "Outcome measures", fill = "Type of Collaboration", 
       caption = "")

ggsave("vign_quota_collab_means_straight.png", plot = last_plot(),
       path = "Figures/", limitsize = F,
       width = 7.5, height = 4.5, units = "in", dpi =300)


## TABLE A.10: Effect of Working with Another Party on Competence Rating: Non-Straightliners Sample, All Treatment Arms --------

# Effect of type of collaboration on competence
vign_h2_comp1 <- lm(vignette_competence ~ vignette_work, 
                    data = mv_straight[mv_straight$vignette_election=="Man/general",])

vign_h2_comp2 <- lm(vignette_competence ~ vignette_work, 
                    data = mv_straight[mv_straight$vignette_election=="Woman/general",])

vign_h2_comp3 <- lm(vignette_competence ~ vignette_work, 
                    data = mv_straight[mv_straight$vignette_election=="Woman/quota",])

# Print table
stargazer(vign_h2_comp1,vign_h2_comp2,vign_h2_comp3,
          title = "Effect of Working With Another Party on Competence Rating: Non-Straightliner Sample, All Treatment Arms",
          font.size = "scriptsize",
          keep.stat = c("n","adj.rsq"),
          model.names = F,
          multicolumn = TRUE)


## TABLE A.11: Effect of Working With Another Party on Pawn-Like Rating: Non-Straightliners Sample, All Treatment Arms ---------

# Effect of type of collaboration on pawnlike
vign_h2_pawn1 <- lm(vignette_pawnlike ~ vignette_work, 
                    data = mv_straight[mv_straight$vignette_election=="Man/general",])
vign_h2_pawn2 <- lm(vignette_pawnlike ~ vignette_work, 
                    data = mv_straight[mv_straight$vignette_election=="Woman/general",])
vign_h2_pawn3 <- lm(vignette_pawnlike ~ vignette_work, 
                    data = mv_straight[mv_straight$vignette_election=="Woman/quota",])

# Output table format
stargazer(vign_h2_pawn1,vign_h2_pawn2,vign_h2_pawn3,
          title = "Effect of Working With Another Party on Pawn-Like Rating: Non-Straightliner Sample, All Treatment Arms",
          font.size = "scriptsize",
          keep.stat = c("n","adj.rsq"),
          model.names = F,
          multicolumn = TRUE)


# TABLE A.12: False Discovery Rate Adjustment Results ------------

## Gather p-values to adjust -----------------

# Use the p-values for the tests of H1.1, H1.2, H2.1, and H2.2
# These are taken from Tables 1 and 2 in the Main Manuscript
# Refer to code in "02_main_manuscript.R"

p_vals <- c(0.966, # H1.1 quota women vs. non-quota women (competence)
            0.583, # H1.1 quota women vs. non-quota men (competence) 
            0.051, # H1.1 quota women vs. non-quota women (pawn-like)
            0.906, # H1.1 quota women vs. non-quota men (pawn-like) 
            0.148, # H2 collab vs. non-collab among quota women (competence)
            0.066 # H2.1 collab vs. non-collab among quota women (pawn-like)
)

## Calculate adjustments and report -------------------------------

bonferroni <- p.adjust(p_vals, 
                       method = "bonferroni", 
                       n = length(p_vals))

holms <- p.adjust(p_vals, 
                  method = "holm", 
                  n = length(p_vals))

bh <- p.adjust(p_vals, 
               method = "BH", 
               n = length(p_vals))

# Table
fdr_tab <- data.frame(cbind(p_vals, bonferroni, holms, bh))

fdr_tab <- fdr_tab %>%
  mutate(
    across(.cols = everything(),
           ~round(.x, digits = 3))
  )

rownames(fdr_tab) <- c(
  "H1.1 Quota Women vs. Non-Quota Women (Competence)",
  "H1.1 Quota Women vs. Non-Quota Men (Competence)",
  "H1.2 Quota Women vs. Non-Quota Women (Pawn-Like)",
  "H1.2 Quota Women vs. Non-Quota Men (Pawn-Like)",
  "H2.1 Co-Partisan vs. Cross-Party Among Quota Women (Competence)", 
  "H2.2 Co-Partisan vs. Cross-Party Among Quota Women (Pawn-Like)"
)

print(xtable(fdr_tab))
